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Abstract 



Using the approach developed in [T] , we succeeded in reconstructing the behaviour of the 
antiferromagnetic Ising model with imaginary magnetic field iO for two and three dimensions 
in the low temperature regime. A mean-field calculation, expected to work well for high 
dimensions, is also carried out, and the mean-field results coincide qualitatively with those 

■ of the two- and three-dimensional Ising model. The mean field analysis reveals also a phase 
structure more complex than the one expected for QCD with a topological 0— term. 

CD ■ 

■ 1 Introduction 



' Quantum field theories with complex actions are systems as interesting as difficult to an- 

J> . alyze: on one hand, the complex action usually gives rise to a severe sign problem, which 

' prevents computer simulations to be performed; on the other hand, the number of known, 

exactly solvable models with complex actions is quite limited. And yet, there are very im- 
portant models in this class, for instance, QCD with a 0— term or at finite baryon density, 
whose solution might lead to an explanation of the strong CP problem and to the under- 
standing of the rich phase structure expected for matter at high pressures. In condensed 
matter physics Haldane showed [2] that chains of quantum half-integer spins with antiferro- 
magnetic interactions are related to the two-dimensional 0(3) nonlinear sigma model with 
a topological term at = tt, and conjectured that such model presents a second order phase 
transition at 6 = tt, keeping its ground state CP symmetric. These are some of the reasons 
^ , why a great effort is being invested in developing new algorithms, capable of dealing with 

complex actions. 

For the particular case of 0— vacuum systems, the partition function in the presence of a 
9 term is periodic and, conveniently normalized, can be decomposed in sectors of different 
topological charge n (or density of topological charge Xn = n/V), as 

Zy (6) = J2Py e''" = E e-^/-(-)e»«^-, (1) 

n n 

which resembles the Fourier transform of the probability distribution function (p.d.f.) of 
the topological charge at 9 = 0. The probability of the topological sector n is therefore 
given by pv (n), and this quantity can, a priori, be measured from numerical simulations. 
Unfortunately, this is a very difficult task, for several reasons: 

i. The precision in a numerical simulation is limited by statistical fluctuations. Thus the 
measurement of pv {n) suffers from errors. 
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ii. Small errors in pv (n) induce huge errors in the determination of Zy {0), as this quan- 
tity is exponentially small, Zy w e^^, due to the sign problem. 

iii. Even if we were able to evaluate pv (n) with infinite accuracy, the terms on the sum 
([T|) differs by many orders of magnitude (from 1 to e~^). 

In fact, the different groups that have tried to determine with high precision the p.d.f. of 
the topological charge either by standard simulations, or by more sophisticated methods (re- 
weighting or multibinning techniques), have found artificial phase transitions in the U{1) and 
CP^ models. The reason behind these ghost transitions is the flattening of the free energy 
for 0— values larger than a certain threshold. In [3J|3], this threshold is roughly evaluated, 
and the flattening behaviour explained. The conclusion is clear: a reliable computation of 
the order parameter (or the free energy) for all values of 9 from the direct measurement of 
the p.d.f. of the topological charge is not feasible due to the huge statistics required. 

That is why other approaches IHl [3 HI IB UHl HH HI] (or at least, serious refinements 
of the standard approach) should be considered. In [T3], a remarkable breakthrough was 
achieved, pushing the threshold of 6 to its limit 6 = tt. The method is based on the 
observation that, since all the coefficients entering in the right hand side of equation ([1]) are 
positive at purely imaginary values of 0, the free energy is given in the thermodynamic limit 
by the saddle point 

/' (x) = h, (2) 

where / (x) represents fv {x„) as 1/ — > oo, and h stand for a purely imaginary 9 field 
9 = -ih. 

Then, the function was fitted to a ratio of polynomials, and integrated analytically to 
obtain f{x). This step is essentially different to what other groups proposed, and it solved 
the problem of the 9 threshold in some systems. Finally, a multi-precision algorithm is used 
to calculate the partition function directly from ([T]), using the function f{x). 

It was demonstrated numerically in IJi, that the errors in the reconstructed / (x) using 
this method were highly correlated. This was a great advantage, for the errors, propagated 
to the exact free energy density, became almost constant, and these errors amounted to 
an irrelevant constant in the free energy. These ideas were successfully tested |13) in the 
one-dimensional Ising model and the two-dimensional U{1) model, and the method was used 
to predict the behaviour of the CP^ model. Furthermore, in reference 'TT the continuum 
9 dependence of CP^ -a confining and asymptotically free quantum field theory- was fully 
reconstructed. Data collapsed for different couplings within percent level and this evidence 
for scaling at non-zero 9 is the strongest indication that the CP symmetry is spontaneously 
broken in the continuum, as predicted by the large N expansion. 

The results were impressive by that time, solving completely the problem of the flattening 
of the order parameter beyond the critical value of 9. The key of this success was the 
aforementioned correlation among the errors: tests performed using the same method and 
adding an apparently negligible 0.1% uncorrelated random error to the measured free energy 
/ (x) led to disaster. Nonetheless, if the error was correlated, it could be as large as 50%, 
and the final result would be quite reasonable. 

However, this method is not yet generally applicable. The flattening can appear -and in 
fact, does- whenever the behaviour of the order parameter is not monotonous. This seems 
to be a general rule. The flattening was first observed in a simple toy model which featured 
symmetry restoration at 9 = n, 

f{9)=\n{l + Acos9). (3) 
For this model, the order parameter vanishes only at ^^ = 0,7r 

A sin 6* . , 

but the method predicted an almost flat behaviour, slightly increasing, beyond the point 
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On the whole, although the method proposed in |13| represented a large improvement 
over what existed at that point, it was clear that another approach was necessary, and that 
is how the method described in the following section was created [Ij . 



2 One-dimensional Ising model 

The Ising model is the simplest model describing ferromagnets, but it is also a good theo- 
retical laboratory to test new algorithms. It is easy to code efficiently on the computer, and 
simulations are very fast, allowing the generation of big statistics, even on large lattices. 
The one-dimensional model in a magnetic field is exactly soluble, which allows us to check 
our results against the exact solution. We can also identify in some sense magnetization 
and topological charge in this model, and regard an imaginary external magnetic field as a 
9 term in the action. Finally, from the numerical point of view, it is even more challenging 
than other complex systems suffering from the sign problem, such as lattice QCD with a 
0— term, yet it remains more accessible. Therefore, it is a good idea to check the goodness 
of any algorithm in this toy model, prior to its application to more physically interesting 
systems. 

The hamiltonian of the one-dimensional Ising model with nearest neighbours coupling J 
and external magnetic field B is 

N N 

H{{si},J,B) = -JY,s^s,+l-BY,s^■ (5) 

i i 

Defining reduced couplings F = J/{KT), h = 2B/{KT)^ the density of free energy is given 

by 

fiF,h) = F + In l^cosh ^ + y e-^P + sinh^ , (6) 

where f{F,h) — l/V\nZ represents the free energy. It is quite remarkable that the Ising 
model within an imaginary external field (i.e., for h — —i9 with S M) is not properly 
defined with ferromagnetic couplings (27j. Setting F > we find that the free energy (O 
becomes undefined for certain values of 0, for the argument of the logarithm may vanish if 
F > 0. Hence, we will deal with the antiferromagnetic Ising model {F < 0) from now on. 

In systems with an even number of spins, the quantity ^ — ^ an integer taking 

any value between —N/2 and N/2. From equation ([6]), the mean density of magnetization 
is 

df sinh| 
(m) = —j- = — . (7) 

^e-^F + sinh^ I 

Equation ([7]) for the magnetization is completely general, in particular it is valid for the 
case of a pure imaginary magnetic field h = i9. For 6 — tt, the Z2 symmetry is restored 
(the magnetic field amounts to a sign a, depending on the 'topological charge' M/2 of the 
configuration a = e"^*^/^; this sign is invariant under a Z2 transformation). Then, the 
question is whether the Z2 symmetry is spontaneously broken or not. Substituting h — > iO 
in dll), we get 

a 

(m) = , (8) 
- sin2 f 

Thus the magnetization takes a non-zero expectation value for the one-dimensional Ising 
model at 9 — TT, a, fact that indicates spontaneous symmetry breaking (see fig. [T|). 



^As the following sections show, the phase diagram of the Ising model within an imaginary magnetic field is 
richer than the one expected for QCD in presence of a ^—vacuum term. 



3 



9 angle 

Figure 1: Magnetization density as a function of the 9 angle for the one-dimensional Ising model. F was 
set to = -0.50. 



The Ising model in one dimension has a striking scaling property. Let us define the 
variable z ~ cosh % , and compute the ratio ^"I'K : 

^ ' ^ tanh -J- 

Therefore y{z) depends on z and F only through the combination (e~^^ — l) ^ z. Due to 
this simplified dependency on F and h, the transformation 



yx{z) = y(^e^z^ (10) 



with A G M is equivalent to a change in the reduced coupling F, or in the temperature 
of the model. The interesting point here is the fact that for negative values of A, this 
transformation can take the variable z = cosh ^ to the range < z < 1, therefore z — cos |, 
corresponding to an imaginary field. This means that we can measure y (z) for imaginary 
values of the magnetic field by mean of numerical simulations at real values of h, which are 
free from the sign problem. 

In order to check if a similar scaling property still holds for other systems, let us assume 
that y(z,F) = y(g(F)z), then 

§^ ^9'iF)z 



1^ 9iP) 



(11) 



To have scaling, the ratio §p/z^ should be independent of h. 

For the one-dimensional Ising model the simulations exhibit a constant ratio over a large 
range of fields h (see Fig. [5]). 

Unfortunately, this property is exclusive of the one-dimensional case. For two dimensions 
the ratio shows a slightly dependence on the reduced magnetic field. For three dimensions, 
the dependence becomes a bit more pronounced. The peak in Fig. [3] is produced by the 
antiferromagnetic-ferromagnetic phase transitior|^. 



^The antiferromagnetic Ising model displays, for strong enough couplings, a phase transition at non-zero 
external magnetic field : the spin-coupling tries to put the system in an antiferromagnetic state, whereas the 
external field tries to order the spins in a ferromagnetic fashion. As the value of the external field increases, the 
ferromagnetic behaviour takes over. 
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Figure 2: Ising's scaling check along formula (fTT|) . The continuous lines represent the analytical result, 
while the crosses stand for the numerical data. We performed short simulations (only ^ 100000 iterations) 
for several values of in a i = 100 lattice. Errors are smaller than symbols. 



In any case, the dependence on the external field for these two models is mild for small 
values of the field h, far from the transition point, and large values of |F| (for low tem- 
peratures). We will see later that this property becomes very relevant when dealing with 
asymptotically free theories. 
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Figure 3: Ising's results in 2D and 3D. The scaling is approximate in 2D and 3D for low values of the 
field. The statistics of the simulations are 100000 iterations, and the lattice lengths are, for 2D L = 50 
and for 3D L = 25. Errors are smaller than symbols. 



3 Computing the order parameter for imaginary mag- 
netic fields 

Although the exact scaling property is absent in higher dimensions, we can still take ad- 
vantage from the methodology it suggests. For the one-dimensional case, a measurement of 
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the order parameter produced at the point {F, z) is equivalent to a measurement done at 
{F' , z') if the fohowing relationship holds 

g{F)z = g{F')z', 

g{F)^{e-^^-iy. (12) 

This way, and choosing carefully the value of a simulation performed at a real value of 
the reduced magnetic field z > 1 is equivalent to another simulation performed at imaginary 
values of h (where z < 1). 

The procedure to find out the order parameter at imaginary values of the reduced mag- 
netic field relies on scaling transformations 1 . We define the function y\ (z) as 

yx{z)^y{e^zy (13) 

For negative values of A, the function y\ (z) allows us to calculate the order parameter 
(tanh ^ y (z)) below the threshold z = 1. If y (z) is non- vanishing for any positive 41, then 
we can plot y\/y against y. Furthermore, in the case that y\/y is a smooth function of y 
close to the origin, then we can rely on a simple extrapolation to y = 0. Of course, a smooth 
behaviour of y\/y can not be taken for granted; however no violations of this rule have been 
found in the exactly solvable models. 

The behaviour of the model dX = tt can be ascertained from this extrapolation. At 
this the model has the same Z2 symmetry as in the absence of field. We define a critical 
exponent 7a 

As z 0, the order parameter tan | y (cos |) behaves as (tt — 6*)^^^^. Therefore, a value of 
7a = 1 implies spontaneous symmetry breaking at 6* = tt. A value between 1 < 7a < 2 signals 
a second order phase transition, and the corresponding susceptibility diverges. Finally, if 
7a = 2, the symmetry is realized (at least for the selected order parameter), there is no 
phase transition and the free energy is analytic at 6* = 7r|f| 

We can take the information contained in the quotient ^ (y), and calculate the order 
parameter for any value of the imaginary reduced magnetic field h = —iO through an iterative 
procedure [1]. The outline of the procedure is the following: 

i. Beginning from a point 'y{zi) — yi, we find the value yi+i such that yx = yt- By 
definition, y,+i = y (^e^ Zi^ 

n. Replace y^ by y^+i, to obtain yi+2 = y (e~^Zi). 

The procedure is repeated until enough values of y are know for z < 1 (see Fig. 0]). This 
method can be used for any model, as long as our assumptions of smoothness and absence 
of singular behaviour are verified during the numerical computations. Indeed the reliability 
of our approach in practical aplications is better when the following two points are well 
realized: 

a. j/(z) takes small values for values of z of order 1. 

b. The dependence on y of the functions yx/y and 7a is soft enough to allow a reliable 
extrapolation. 

In the one-dimensional model these two properties are realized in the low temperature regime 
(see equation ini and Fig. 5), but the two and three-dimensional models, notwithstanding 
they do not verify a perfect scaling law as in the one-dimensional case, they also show a 

^Even though the possibility of a vanishing y (z) for some value z > can not be excluded completely, it does 
not happen for any of the analitically solvable models we know. 

^Other possibilities are allowed, for instance, any 7a > 1, 7a G N leads to symmetry realization for the order 
parameter at 6 = tv and to an analytic free energy. If 7a lies between two natural numbers, p < < q, p, q £ N, 
then a transition of order q takes place. 



6 



very good behaviour (see Figs. 7, 8). Indeed the relevant feature, at least in what concerns 
point a, is that, at low temperatures, the magnetic susceptibility at small values of the 
real external magnetic field h, which is essentially y{z), takes small values; and this is also 
true for any dimension. In the more interesting case of asymptotically free models, the 
analogue of the magnetic susceptibility is the topological susceptibility, and it is well known 
that topological structures are very suppresed near the continuum limit. Therefore, and 
on qualitative grounds, we expect a much better implementation of our method in the low 
temperature regime of the Ising model and near the continuum limit of asymptotically free 
theories. A check of this statement for the Ising model is the content of this article, and 
concerning asymptotically free models, the method was successfuly applied to the analysis of 
the continuum 0— dependence of CP^ |14j . showing a very good realization of points a and 
b. In the more general cases we should find out whether the model complies with these two 
points or not, and pleasant surprises are not excluded. Indeed, we checked in [28 that the 
conditions for the aplication of our approach to CP^ also hold, and this allowed us to verify 
the Haldane conjecture and the relevant universality class of the non-linear 0(3) tr-model 
in two dimensions. 




0.00 0.05 0.10 0.15 0.20 0.25 

y 

Figure 4: Iterative method used to compute the different values of y{z). y\ is ploted as a function of y 
using a dashed line in the region where direct measurements are available, and a continuous line in the 
extrapolated region. The straight continuous line represents y\ — y. 



4 Numerical results 

The first thing we did was to test the method in the one-dimensional Ising model, and we 
checked the results against ([S]). The simulations were performed at a fixed volume, N = 1000 
spins, and fixed reduced coupling F — —2.0. As the one-dimensional Ising model has no 
phase transitions, and furthermore enjoys the exact scaling property, there is no point in 
checking the method for several values of the reduced coupling. The parameter we varied was 
the reduced magnetic field h. As the simulations were quite fast, we could obtain data for 
many values of h with large statistics. In fact, for each point in the plots we performed lO'' 
metropolis iterations. In order to reduce autocorrelations, we performed at each iteration 
two sweeps over the lattice, proposing metropolis changes in the spins. The plots for the 
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Figure 5: Calculation of the critical exponent 7a. The crosses correspond to the numerical simulation 
data, whereas the line is a quadratic fit. The extrapolation to zero seems quite reliable, as the function 
is smooth enough. Errors are smaller than symbols. 
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Figure 6: Order parameter as a function of 9. 
spontaneous breaking of the Z2 symmetry at 9 



The non-zero value of the order parameter marks the 

= TT. 



critical exponent and the order parameter are shown in Figs. [S] and [51 Our result for the 
critical exponent from the fit in fig. [S]is 

7A = 0.99980 ±0.00008, 

which agrees with the analytical result. 

Then we simulated higher dimensional models, expecting to see departures from this 
behaviour, as these models feature phase transitions between ordered (antiferromagnetic) 
and disordered phases. The two-dimensional simulations were done in a 100^ lattice, after 
100000 termalization sweeps. We spent 5000000 steps to measure each point accurately. 
The three-dimensional case, on the other hand, used a 50"^ volume, and measured each point 
for 2500000 steps after 100000 steps of thermalization. The results showed the expected 
departure in the behaviour. Our result for the critical exponent 7^ w 2 reveals a vanishing 
order parameter at 6' = tt in the ordered phase {F — —1.50 for 2D and F = —1.00 for 3D), 
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as shown in Figs. [71 HI 



= 1.9997 ±0.0002 
= 1.9998 ±0.0002 



We can confirm this facts by plotting the order parameter against 
andfni 



2.10 



as it is done in Fig. [10] 
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Figure 7: Calculation of the critical exponent 7^ in the ordered phase of the two-dimensional model. 
The pluses correspond to the numerical simulation data, whereas the Hne is a cuadratic fit. Errors are 
much smaller than symbols, except for the points lying close to the origin. 



2.4 



2.3 
2.2 
2.1 
fl 2.0 
1.9 



1.8 - 
1.7 
1.6 



3D Ising Model F = -1.0 i-H — 1 
Constant Fit 



»)H) H ) U4 tlll> - > - III lllli mn M I I I — HH 1- 



5e-05 



0.0001 



0.00015 



0.0002 



0.00025 



y 



Figure 8: Calculation of the critical exponent 7a in the ordered phase of the three-dimensional model. 
The pluses correspond to the numerical simulation data, whereas the line is a constant fit. Errors are 
much smaller than symbols, except for the points lying close to the origin. 



The disordered phase revealed a caveat of this method, as it was impossible for us to 
extrapolate the function ^ (y) to zero. The reason is simple: at small values of y 
and yx approach unity, for at vanishing F we recover the paramagnetic Langevin solution 

^Actually the Z2 symmetry is spontaneously broken, for the staggered magnetization ms 7^ [25]. This point 
will be clarified in the mean-field approximation. 
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Figure 9: Due to the peaked behaviour of the quotient y\/y around y — 0.3, the extrapolation to zero 
required many simulations at small values of the magnetic field. Errors are much smaller than symbols, 
except for the points lying close to the origin. 
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Figure 10: Order parameter as a function of 9 in the ordered phase of the two-dimensional model. It 
vanishes at 9 — n. 



m — tan I and y = 1. The smaller the value of F, the greater the gap between zero and 
our data becomes, and at some point, the extrapolation is not reliable any more, and the 
results depend strongly on the fitting function used. An example can be seen in Fig. 1121 
where the data for the two-dimensional model at F = —0.40 are plotted. In this case, we 
are too far from zero to find out accurately the critical exponent, and the value of F could 
not be lowered much more, for the transition to the ordered phase is known to happen at 
F —0.44. In Fig. [T3] a similar example is shown for the ordered phase in the three- 
dimensional model, but this time a tentative extrapolation could be done, casting a reliable 
result. 

These examples show how this method works fine when the antiferromagnetic couplings 
are strong enough. In general, as discussed at the end of the previous section, the method 
performs well in the low temperature regime od the Ising model and for asymptotically free 
theories, whose continuum limit lie in the region of weak coupling. In this region, the density 
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Figure 11: Order parameter as a function of 6 in the ordered phase of the three-dimensional model. As 
in its two-dimensional counterpart, it vanishes at 6 = tt. Errors are smaller than symbols. 
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Figure 12: Failed calculation of the critical exponent 7a in the disordered phase for the two dimensional 
model. Our data is so far from the y = axis, that an extrapolation can not be used to find out the 
value of 7a. A peak for lower values of y, as the one appearing in Fig. [HI cannot be discarded 'a priori'. 
Errors are much smaller than symbols. 



of topological structures is strongly suppressed. Thus the order parameter and y (z) take 
small values, making the plot ^ (y) easily extrapolable to zero. In the particular case of the 
antiferromagnetic Ising model, large values of \F\ ensure a small magnetic susceptibility. A 
high value of the dimension also helps, for instance, the three-dimensional model requires a 
lower value of the coupling than the two-dimensional case to make a reliable extrapolation 
of ^ (y) to y — 0, for each spin is affected by a higher number of neighbours. 

As this method failed to deliver interesting results in the disordered phase, we tried 
a different approach: we expected naively that the two-dimensional model resemble the 
one-dimensional model at low values of the coupling. Since the p.d.f. method worked well 
for the one-dimensional case |13) . it made sense that we applied it to the present scenario. 
What we found is an unstable behaviour: sometimes the method seems to predict the phase 
transition, in the sense that at finite volume there is not true phase transition, and an abrupt 
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Figure 13: Another calculation of the critical exponent 7^ in the ordered phase F — —0.3 for the three- 
dimensional model. Our data approaches the y = Q axis enough to try an extrapolation, but the result 
suffers from much larger errors than in the F — —1.0 case. Here 7a = 2.079±G.003, and the measurement 
errors are much smaller than symbols. 



modification in the order parameter, linking the two expected behaviours, should happen. 
This is what we observe in one of the data sets of Fig. [141 Nonetheless, if a slightly different 
set of points is taken to fit the saddle point equation ^ , the resulting data show a sharp 
departure from the expected behaviour at some 9. 
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Figure 14: Failed calculation of the order parameter in the disordered phase using the improved p.d.f. 
method described at the beginning of this paper. In the first case, the points seem to predict a phase 
transition at 9 ^ 2.35, whereas in the second case the points depart sharply from a smooth function 
and never come back. The only difference between fits was the number of points used: in the second 
case, only half of the points (the closest to the origin) were used. Other modifications in the fitting 
procedures indicate us that the transition point is not stable. This might indicate either a failure in 
the fitting function, or a phase transition, and the impossibility for the method to precise the transition 
point, unless a perfect ansatz is made. Errors were not estimated. 



There are two possible explanations to this behaviour: either the fitting function selected 
is wrong, or there is some hidden phenomenon we are overlooking. The fitting function used 
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was an odd quotient of polynomials 

ax^ — X 
cx'^ — b 

which should account well for the behaviour of the order parameter, given the assumption 
that it is similar to the one-dimensional case. The addition of more terms to the fit did not 
do much to improve the result, hence this possibility was discarded. 

The existence of a phase transition in the middle, however, was an interesting option. 
Indeed, the two-dimensional model in the presence of a term was solved exactly at the 
point 9 = n almost sixty years ago by Yang and Lee in [SS], and reviewed again in [50] . 
In those papers, a phase diagram was proposed were the antiferromagnetic model always 
stayed in an ordered phase at any non-zero value of F. Since the system is in a disordered 
state for low _F's and zero field, some phase transition has to occur in the middle. Thence, 
the failure of the p.d.f. method should be due to a poor ansatz for the fitting function, 
caused by the presence of a phase transition at some 9c. 

The fact that the results for the two- and three-dimensional models are qualitatively 
the same in the ordered phase, makes us wonder whether this behaviour changes for some 
value of the dimension D > 3. Moreover, the behaviour of this model in the disordered 
phase is unknown to us. That is why we decided to carry out a mean-field approximation of 
the model, and compute the critical exponent 7a. As we know, mean- field results for other 
critical exponents are exact for the n-dimensional ferromagnetic Ising model, provided that 
n > 4. Thus we expect that, if the mean-field result for 7^ is the same to that of the two- 
and the three-dimensional Ising model, then "fx — 2 for any value of the dimension. 



5 Mean-field calculation 

In antiferromagnetic compounds the spin- alignment pattern is staggered, hence, in order 
to define the mean-field version of the antiferromagnetic Ising model, we should divide the 
lattice in two sublattices, and define a coupling among spins whose sign depends on whether 
these two spins are on the same sublattice or not. For spins belonging to the same lattice, the 
coupling should be ferromagnetic (J > 0), but for spins belonging to different sublattices, the 
coupling should favour antiparallel ordering (—J < 0), according to the antiferromagnetic 
nature of the system. Therefore, two different mean-fields should appear, (si)igs^ — nii and 
{si)ieS2 = '^2, referring to each one of the different sublattices. 

A reasonable mean field hamiltonian compatible with these requirements is H 




H(J,S,{sJ) = -- > s.-> -Bi > s,-i?2 > s,. (15) 



We define as before hi = (for each sublattice i — 1,2) and F — and the usual 
and the staggered magnetizations, m = mi -I- 7712, ms — mi — m2- We are interested in the 
model for imaginary values of the field, h — > i9. Using standard saddle point techniques (see 
appendix A for details) we obtain the mean field equations: 



2 cosh2(2|_F|ms)-sin2 | 

1 sinh(4|_F|mg) 

2 cosh^(2|F|m5)— sin' 



TO^ = i sinh(MF\ms) QyX 



An analysis of these equations shows that for low values of F there's only one solution, 
7715 = 0, corresponding to a paramagnetic phase. For values over a certain Fc the situation 



®Of course other mean-field approaches could be used, but we expect to obtain the same qualitative picture 
for the phase diagram. For example, a calculation in the scheme of [31] produces similar results, and more 
importantly the same value for 7a as the one in our calculation. 
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changes and two new symmetric solutions appear, which are in fact the physically relevant 
ones. From the saddle point equations, pop and ([5T|). Fc can be obtained as a function of 6: 



2Fr = cos"' 



(18) 



The phase diagram of the system in the F-6 plane is shown in Fig. [151 There is a second 
0.60 



0.50 




0.40 - 



fcn 0.30 - 



0.20 - 



0.10 - 



0.00 



Figure 15: Phase diagram of the mean- field approximation to the antiferromagnetic Ising model in the 
F — 9 plane. 



order phase transition at the critical line (1181) . As — > tt the paramagnetic phase narrows, 
until it is reduced to the single point F = aX = tt. The staggered phase (with antiparallel 

cos^ - 

spin ordering ms ^ and F > — 2^) 011 the other hand features Z2 spontaneous symmetry 
breaking at = tt, as equation [T7] shows. The fact that this model features a phase tran- 
sition at non-zero values of the external field is quite remarkable. This kind of transitions 
would never appear in a ferromagnetic model within a real external magnetic field, as the 
external field and the spin coupling work in the same direction: parallel spin alignment. On 
the contrary, in the antiferromagnetic case, the introduction of a real external field produces 
frustration, whose origin comes from the competition of the spin coupling, trying to move 
the spins towards an antiparallel configuration, and the external field, favouring a completely 
parallel structure. What is remarkable of these results is the fact that in the antiferromag- 
netic case an imaginary magnetic field, strong enough, is able to move the system from the 
paramagnetic phase to a phase with antipararell ordering. 

All the magnetizations are continuous functions, but the staggered susceptibility xs di- 
verges, as usually happens in a second order phase transition. The topological^^ susceptibility 
XT, on the other hand, displays a gap at the critical line. The computation (see Appendix 
A) gives the result 



lim XT 



lim X.T 



3 2|F|-1 
4|F| 4|F| - 3' 



(19) 



Finally, the critical exponent 7^ for this mean-field theory can be calculated, to see 
if it coincides with that obtained in simulations. In order to do so, we expand m in the 
neighbourhood oi 9 — tt 



m{6) ^ m (tt) 



dm 
~d6 



(tt - 0) -f 



ff^m 



(9612 



[tt - ef + 



(20) 



^Topological in the sense that M/2 is a quantized charge, m/2 is its associated charge density, and xt the 
susceptibility. 
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If 7a is not natural number, we expect the first non-zero derivative to diverge. On the 
contrary, if 7^, is a natural number, the order of the first non- vanishing derivative will give 
us the critical exponent. Taking derivatives 



dm 
'do 



(2 cos^ e~l) cosh^ (2 \F\ ms) + I - cos' 



2 £ 
2 



2 (cosh2(2|F|ms)-sin2 
i2|f|singsinh(4|F|ms) ^|,^, 



dms 



(cosh^ (2 |F| ms) 



,2 e^ 



sinh^ (2 |F| ms) 



+ 



2|F|sin0sinh(4|^^|ms) ^|,^ 



sinh* (2 |F| TOs) 



(21) 



The first term on the r.h.s. of ^ does not diverge since ms is not vanishing as 9 ^ tt. The 
second term is proportional to 

hm smw. 

After a tedious calculation, it can be shown that it vanishes, therefore 



m{d) 



2sinh^ {2\F\ms) 



(22) 



with K a non-zero constant. The magnetization behaves as (tt — 9)'^'"~^ in the neighbour- 
hood of 6' = TT. Thence, for mean-field antiferromagnetic theory 7^ = 2 for F ^ 0. For 
^ = 7a = and the magnetization diverges as tan | when 9 ^ n. 

Since mean-field theory works better in high dimensional systems (it reproduces all the 
critical exponents exactly for the ferromagnetic Ising model in dimension 4 and above), and 
the exponent 7a seems to have settled in 7a = 2 for the two- and three-dimensional models, 
and for the mean-field approximation, we expect this result to hold for any dimension of the 
system. This is not a proof, but in fact, it would be very remarkable if the behaviour of the 
antiferromagnetic Ising model in a higher dimension departed from 7a = 2. 



6 Conclusions 

Although the aim behind this investigation of the antiferromagnetic Ising model is to test 
our techniques for a future application to QCD in presence of a 6* term, the results obtained 
through this work deserve attention on their own merit. Using the method described in 
section[3l the order parameter for the Z2 symmetry can be calculated for any value of 9, and 
although there are some regions of the phase diagram where the method does not work well 
(high temperature regime), it provided us with enough information to make an educated 
guess on the phase diagram of the theory. 

Our guess was later confirmed by a mean-field calculation, which shares many properties 
with the original model. The results of [29] and [30] supplied the remaining information for 
the two-dimensional case. In the end, we were able to reconstruct qualitatively the whole 
phase diagram of the theory for two-dimensions, and although we did not pursue to solve the 
model for higher dimensions, the mean-field calculations, and the fact that the behaviour 
for the two- and the three-dimensional models is the same for large values of F, give us 
strong indications that this phase diagram is qualitatively valid for any dimension of the 
model larger than one. The one-dimensional model is an exceptional case in which only 
one phase appears with spontaneous magnetization at 6 = tt. On the contrary for d = 2, 3 
our numerical simulations show a density of magnetization that continuously vanishes at 
9 = TT at low temperatures, and the mean field calculation strongly suggests that this result 
holds for any temperature and any larger dimension. However this does not mean that 
the Z2 symmetry of the model at = tt is realized in the ground state since the mean 
field calculation shows a non vanishing value of the staggered magnetization at 9 = tt. 
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Indeed there is a region of non- vanishing measure in the F — 9 plane, including the 9 = ir 
line, where the staggered magnetization does not vanish and the saddle point equation (|17p 
shows two symmetric solutions for this quantity. In all this phase translational invariance is 
spontaneously broken and in the 9 — tt line parity is also spontaneously broken. 
The method only has two caveats: 

i. It does not work properly (and can give wrong results) if there is a phase transition 
for some 9 < tt. 

ii. For small absolute values of the coupling F (high temperatures) the required extrap- 
olations are not feasible. 

Fortunately, the standard wisdom on QCD, based on reasonable assumptions, expects 
only one phase transitions at 9 — tt [32], and QCD is an asymptotically free theory, thus its 
continuum limit lies, as discussed in section 3, in the "low temperature" regime where our 
approach works well. Therefore this method has become the perfect candidate to perform 
simulations of QCD with a 9 term, which might provide precious information concerning 
the strong CP problem. 

A final remark on the antiferromagnetic Ising model at = tt is pertinent. By using 
polymerization techniques an algorithm able to simulate the model at 9 = tt and free from 
the sign problem can be developed, and it could be useful to test if the mean field predictions 
reported here are verified for any dimension larger than 1 and any temperature at 9 = tt. 
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A Mean field model 

The hamiltonian of the mean- field model is: 

H (J, B, {s.}) = f E - E *0 " ^1 II - ^2 E 

\ieSi jeS2 / leSi jeS2 



(23) 



where we have assummed that an independent field Bi acts on each sublattice i — 1,2. This 
modification will allow us to compute separately mi and TO2 once we have obtained the free 
energy; but of course, for the case of an uniform external field we shuld set B = Bi = B2 
at the end of the calculation. 

We define as before hi = i = 1, 2 and F — -j^, and the usual and the staggered 

magnetizations, m = mi -I- TO2, "ms ~ fni ~ ^2- The partition function 



Z{F,h)^Y.^ 



-hi. E 



h2. E 



2 jeS2 ' 



(24) 



can be summed up by applying the Hubbard-Stratonovich identity to linearize the exponent 

(25) 



Z{F,h)^^l E^ 













2x^-h2 


N 2 




N 2 



jeS2 "3 



dx. 



At this point we see that the introduction of the 9 term through the transformations 

hi -J> i9i, h2 i92, 
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render the hyperbolic cosines complex. The ^ factor allows us to define properly the quan- 
tized number The integrand factorizes, as there is no spin-spin interaction 



■)N r°° 



Z{F,h) = — 
7r2 



cosh 2x- 



\F\l 



cosh 2x — i — h i — 
\ N-2 2 , 



dx. 



(26) 



Now we bring the transformation 



X N 2y 

dx N^dy 



so (|^ becomes 



Z{F,h) 



i \n\cosli[2\F\2 y+i^) cosh( 2 1 F| 2 i 



N 



dy. 



(27) 



where we have written the whole integral as an exponential. We can evaluate the integral 
in the large N limit using the saddle-point technique. 

1 



-|- lim — In 

Af->oo N 



lim — In Z (.LB) = ln2- 



-y +^ In cosh^(2|F|2 j/)-s 



N 



dy. 



The maximum of 



gives us the saddle-point equations 



yo 



■In 



cosh^ (2\F\'^ y'^ - sin^ 



sinh(4|F|5yo 



cosh^ 2 |F| ^ yo - sin 



-1 + 2|F| 



Thus, the free energy is 



cos^lcosh (a\F\'^ 



yo 



sinh' 



(2|Fp 



yo 



cosh^ (2\F\'^ yo) - sin^ f 
/(F,/i)-ln2 + g(yo) 



0, 



< 0. 



where yo verifies (150]) . (|?T|) . 

We can also evaluate the magnetizations: 



mi 



7712 = 



77T. : 



^ cosh(^2|F| 2 sinh(^2|F| 2 j/o)+j sin f cos f 
2 cosh2(2|F|3yo^-sin2 I 

cosh(2|F|5yo) sinh(2|F|5j;(,^-isin | cos J 
^ cosh='(2|F|5)/o)-sin2 | 

i sin ^ cos ^ 



777s 



cosli2(2|F|7y„j_sin2 I 

cosh(2|F|^i/o) sinh(2|F|3yo) 
cosh^(2|F|ii;o)-sin2 | 



(28) 
(29) 

(30) 
(31) 

(32) 

(33) 

(34) 
(35) 

(36) 
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Therefore, and using (PO)) . 

yo^\F\^{ms). (37) 

In order to compute the gap in the susceptibihties at the critical hne, we must examine 
the behaviour of the variable j/o in the neighbourhood of 9c- The way to proceed is to expand 
the hyperbolic functions in ((30|) as a power series in yo. For 9 < 9c the only solution to the 
saddle point equation is j/o = Oj so we expand around this point 

smh[4\F\^yo) = 4 yo + y ^o' + O (y^) , (38) 
cosh (2 |Fp 2/0) = 1 + 2 |F| yl + O (y^) , (39) 
so the saddle-point equation becomes 

4 |F| y^ + cos^ I 
Now we expand again the denominator of (j40p up to j/g, 

,1^1 2^, — ^ = - -^^yo + o (2/0) • (41) 

4 |F| J/o + cos^ I cos^ I cos^ I 

Therefore, for 9 > 9c, 

2 \F\ , ,2 2cos2 I - 3 . 

2/0 ^ ^|yo + 8 ^ (42) 

cos^ I 3 cos* I 

We already know of the yo = solution. Solving the quadratic equation that is left. 



3 (cos2 f - 2 



cos 



2 



Thus yo tends to zero as 9 approaches the critical value for a given F . Its derivative with 
respect to 0, on the other hand, 



dyo _ / 3 sin 9 (cos4 f - 3 cos^ f + 3 |F|) 



diverges as 



2 cos2 f - 3) [(cos2 f - 2 |F|) cos^ |] 



(44) 



2^-2 '^"'^2 . 



cos^ ^ - I \t \] cos 



2 ^ 1^ ly ^"■^ 2 

for at the critical line cos ^ = 2F. The divergence cancels in the product yo^- As 
yo — \fW\'^s, this also applies to Tig^^. 

The solution obtained in (03]) can be used to calculate the behaviour of the susceptibilities 
around the critical point. The 'topological' susceptibility 

dm dm d9 (cos2 | - l) cosh (2 |F| ms) + 1 - cos2 f 

XT 

'"''2 ""'2 



di% d9 dil, {cosh^ {2\F\ms {9)) 



2|F|sin6lsinh (4|F|tos)^ 
(cosh^ (2 |F| ms) — sin' 



2 B 
2) 



,2 ' 



(45) 



takes the value 



lim XT = — \-f- = TT^, (46) 

9^67 C0S2 ^ 2 li* I 
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as we approach 9c from below. However, if we come from the antiferromagnetic phase 9 > 9c, 
the second term gives a non-zero contribution, for the derivative diverges at the critical 
line. The divergence is cancelled exactly by the factor sinh (4 \F\ ms), as explained before, 
and what remains is a finite contribution 



dm,'- 



ms- 



d9 



3 sin 9c 



16|Fr(4|F| -3) 



2 I I sin 6* sinh (4 |F| mg) ^ 
(cosh^ (2|i^|TOs) -sin^ f)^ 

1 



In the end 



and the gap is 



3 sin 9c 



lim yt , , 



32|Fr(4F-3) 
3 2|i^|-l 



41F| 4|F| - 3 

3 2!F|-1 



(47) 



(48) 



The staggered susceptibility diverges at the critical line. This is quite expected, as for 
9 = the susceptibility diverges at the critical point. In order to obtain xs, we need to take 
derivatives with respect to a staggered field 9s, and then take the 6*5 —J> limit. To this 
purpose, we use the original form of the free energy (f^ with a 9s term 



f{F,ms,9,9s) = -\F\ml + ^\n 



cosh 2 \F\ms 



cosh 1^2 |F| ms 

i9-ies^ 



i9 + i9s 



(49) 

Taking derivatives with respect to mg we should recover the saddle-point equation with the 
addition of the ^5 source 



df 
dms 







-2\F\ms + \F\ 



tanh 2 \F\ ms 



i9 + i9s 



+ 



tanh 2 |F| ms - 



iW - it's 



(50) 



From (f5| a new equation for the staggered magnetization is obtained 

i9 + i9s'' 



ms 



tanh 2 \F\ ms 



tanh 1^2 |F| ms - 
The derivative with respect to 9s gives us the susceptibility 

dms 1 + 2 |F| xs 



(51) 



XS 



dl^ 2cosh2(2|F|ms + ^^^^) 



1 + 2|F|X5 



1 + 2|^^|XS 



where 



X = 



2cosh^ {2\F\ms 
1 



1 



X, 



(52) 



2cosh^(2|F|m5 + i^^) 2 cosh^ (2 |F| 



Moving all the terms proportional to xs lo the l.h.s. 

2xs(l-|i^|^)=^, 



(53) 
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we can find the value of ^5 



X 



Xs 



2-2 F X 



(54) 



Tlie quantity X must be evaluated at the point 9 — 9c and 9s — This is not a difficult 
task and the final value is 

The staggered susceptibility, on the other hand, diverges at the critical line. Approaching 
the critical point from the high temperature region (where the only solution to the saddle 
point equation is ms =0), we find 



1 



XS 



2\Fr\ -2\FV 



(55) 



the susceptibility diverges at the critical line F = Fc and the critical exponent for this 
divergence 

xs«|r-r,r^^ (56) 

is 7s = 1. 

Finally, and to elucidate the behaviour of m (9) as 6 ^ tt, we need to work out the 
following limit 



lim 



d9 



■ sm ( 



As sin 9 — > when 9 approaches tt, only if the derivative diverges at 6* = tt is the 
product (j57p non- vanishing. The expansion we performed previously is not very useful here, 
as the point 6* = tt is far from the critical line (unless we are taking the F ^ limit as well). 
The way to solve this problem is to compute implicitly the derivative from the saddle-point 
equation (IT7|) at = tt 



dm 5 



d9 



dms 2|F|cosh(4|F|m5) 



d9 cosh2(2|F|ms) - sin^ 
sinh(4|F|ms) ^ ^ , _ . 



^ |F| sinh (4 |f I ms) 



sin 9 

4 



dms 



(cosh^ (2 |F| ms) -sin^f)^ 

4 \F\ cotanh (2 |F| ms) [1 - cotanh (2 \F\ ms) 
Moving all the terms to the l.h.s. of the equation we find that either 

1 - 4 |F| cotanh (2 \F\ ms) [1 - cotanh (2 \F\ ms)] = 0, 



d9 



or 



dms 



d9 



0. 



(58) 



(59) 



(60) 



The first case is impossible, for the solution to the saddle-point equation at 9 = tt imposes 

™s (t^) = cotanh (2 \F\ ms) , (61) 

which is non-zero and verifies 

\ms\ > 1, 

so the l.h.s. never vanishes, for the second summand is always positive. Therefore, 
applies and the derivative vanishes at = tt. 
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